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We propose a simple and general model accounting for the dependence of the viscosity of 
a hard sphere suspension at arbitrary volume fractions. The model constitutes a continuum- 
medium description based on a recursive-differential method that assumes a hierarchy of 
relaxation times. Geometrical information of the system is introduced through an effective 
volume fraction that approaches the usual filling fraction at low concentrations and becomes 
one at maximum packing. The agreement of our expression for the viscosity with experiments 
at low- and high-shear rates and in the high-frequency limit is remarkable for all volume 
fractions. 

PACS numbers: 



I. INTRODUCTION 



A large amount of work has been devoted to the understanding of the rheology of suspensions 
due to its ubiquitousness and the central role that they play in many technological processes. To 
control the structure and flow properties of such suspensions is of crucial practical importance. 
Due to its many body character, the study of the rheology of these systems is extremely difficult, 
even the simplest suspension, composed of identical hard spheres, show rather complex rheological 
phenomena. Since Einstein's calculation of the viscosity for a dilute suspension of spheres [1[ , a lot 
of effort has been devoted to the extension of Einstein's result to more concentrated suspensions in 
which correlations among particles must be taken into account. However, in spite of the remark- 
able works which make central contributions to the understanding of the rheology of concentrated 
suspensions Qj-lll], the problem still constitutes a challenge for theoretical description [16]. Es- 
sentially, most of these works introduce particle correlations by taking into account hydrodynamic 
interactions, and provide a conceptual framework that explains how these interactions modify the 
behavior of the viscosity as a function of the volume fraction <p for arbitrary values of the wave 



2 



vector and the frequency: rj(k,uj;(f)). Sometimes, these approaches provide compact expressions 
for the effective viscosity, although usually a multipole expansion must be done in order to obtain 



m. 



numerical values 

Beenakker {J, by neglecting Brownian motion, gave a theory for tin* concentration and wave- 
vector dependence of the viscosity valid up to high concentrations that fully takes into account 
the many-body hydrodynamic interactions between an arbitrary number of spheres. Batchelor and 
Green P] developed a theory for the evaluation of the bulk stress in a semi-dilute suspensions of 
hard spheres. By neglecting Brownian motion and assuming a random particle distribution they 
found the following expression for the viscosity 

v (<P) = m U + \<\> + 5.2^ , (1) 

where r/o is the solvent viscosity and (j) is the volume fraction. This result has been regarded as 
the high-frequency value of the real part of the complex shear viscosity 12| . A different result is 
obtained if Brownian motion is taken into account. Then appears an extra term accounting for 
the stresses generated in the dispersion by the random movements of the particles. Numerical 
calculations performed in Ref. || lead to the formula 

r / (<A) = %fl + ^ + 6.170 2 V (2) 
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considered as the low-shear and zero- frequency value of the real part of the complex shear viscosity 



12] 



because in the calculation Brownian motion dominates 

Bedeaux and collaborators [7|] introduced a general formalism for the calculation of the effective 
shear viscosity for a suspension of spheres as a function of the wave vector and the frequency 
for arbitrary values of the volume fraction. Using the conservation law for the total momentum 
density of the fluid and the suspended spheres, they perform a multipole expansion of the one sphere 
response function to second order in the wave-vector and an expansion in the volume fraction up 
to the second order, obtaining the expression 

V (<t>) =VoU + ~<l> + 4.8^ 2 ) . (3) 

Russel and Gast [10] formulated a theory for the nonequilibrium structure and stresses in a sheared 
suspension with a fluid rest state. Many body interactions were handled exactly in the thermody- 
namics but truncated at the pair level for the hydrodynamics. 

More recently, on the basis of the pair Smoluchowski equation, Cichocki and Felderhof 29| have 



obtained exact results for rj(cf>,uj) to 0(4> 2 ). Without Brownian motion contributions, their result 
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for u = is 

7/(0) = fjofl + -0 + 5.00^, (4) 

while with Brownian motion contributions they find 

V(4>) = VQ U + ^ + 5.910 2 V (5) 

Verberg and collaborators [2] presented an approximate theory based on the contribution of two 
physically relevant processes taking place at different time scales: (i) at short time scales the 
viscosity of the suspension increases when compared to that of the pure solvent at infinite dilution 
due to the finite probability to find two particles at contact, whereas (ii) at long times it is difficult 
for a Brownian particle to diffuse out of the cage formed around it by its neighbors. For uj = 
these authors obtain the approximate result 



while for uj — > 00 they propose 



r?o(l + ^ + 6.03</» 2 ) , i ) 



r r\ 1-0.50 

,M = m ^ w , (7) 



which is accurate for all < < 0.55. Although their results compare reasonably well with exper- 
imental data, they are based in a number of assumptions not completely justified, as mentioned in 
Ref. 

The simple formulas reviewed are of fundamental importance for the understanding of the 
rheology of suspensions. However, these formulas are far from satisfactory since they only lead to 
good quantitative results in the low volume fraction regime 



This restriction have lead to many works searching for simple procedures or semi-empirical ex- 
pressions to calculate the viscosity of suspensions and dispersions at frequencies and shear rates 
relevant to industrial applications [l^]. Several phenomenological formulas [3] have been pro- 
posed in order to fit experiments in the largest possible range of volume fractions. For example, 
Krieger and Dougherty found empirically [18f] that the expression 

V^) = Vo(l-4-) (8) 



with <j) c the filling fraction at maximum packing, agrees reasonably well with the experimental data. 
Moreover, this relation reduces to the correct Einstein equation in the limit of infinite dilution. 



Other example is the following expression obtained by Quemada [1 

r) (<!>)= ^f 1 -^) > ( 9 ) 

which gives the right divergence for 77 as approaches </> c , but, it does not reduces to Einstein's 
expression for low (j>. Nonetheless, although useful, these expressions lack of physical support and 
therefore do not provide insight about the phenomenology of the system. 

In this article, we explore a different approach to the problem leading to a simple an quantita- 
tively powerful description of concentrated suspensions of mono- and poly-disperse suspensions of 
colloidal particles. Our model assumes that at macroscopic length and time scales, the suspension 
can be considered as a truly continuum medium. However, the precise formulas are derived only 
by considering the discrete nature of the suspended phase through its possible configurations and 
poly-dispersion. 

The technique we will adopt has been introduced long time ago in the context of different 



physical properties of dispersed systems 



2i| 



23, 



24 



25[ and nonlocal dielec- 



22l | . including optical 

trie systems [26]. In the case of rheological properties, a limited success has been obtained when 
comparing wit h exper imental results because the geometrical information of the system has poorly 

28^ . In this work, we improve previous results by incorporating struc- 



been included 
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tural information from the start through an effective filling fraction. Moreover, in order to better 
understand some of the assumptions implicit in the differential technique, we perform a mesoscopic 
analysis of the problem by using Fokker-Planck equations 



31 



32 



34( | . This analysis allows us to 



elucidate the main assumptions and approximations underlying the differential method. 

The article is organized as follows, in Section II we introduce our differential effective medium 
model. In Section III we perform an analysis of the model based on a mesoscopic and hydrodynamic 
approach. In Section IV we compare the predictions of our model with previous theories and with 
various experimental results. Finally Section IV is devoted to conclusions. 



II. DIFFERENTIAL VISCOSITY MODEL 



At very low volume fractions, the shear viscosity r] of a suspension of hard spheres can be 
predicted by the simple formula [l( 

v(<t>) = VoU + lA, (10) 

that Einstein calculated from the viscous dissipation produced by the flow around a single sphere. 
Notice that expression (|10p does not take into consideration the way in which the spheres are 
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arranged in the sample. Therefore it contains no information about the maximum packing of 
spheres the system may allocate. This information is not relevant for very low volume fractions, 
nonetheless it does matter for higher concentration of spheres. Even more, Einstein's expression 
contains no information about the radii of the spheres and it may be applied to any polydisperse 
distribution of spheres at low volume fractions giving to it its universal character. 

In order to extend Einstein's expression to larger filling fractions <p we use a differential effective 
medium approach. For this purpose, it is important to consider that experimental results show 
that the viscosity diverges at volume fractions (j) < 1. This is related to the fact that the spheres 
can not occupy all the volume of the sample due to geometrical restrictions. For example, for 
a face centered cubic (FCC) arrangement of identical spheres, the maximum volume that the 
spheres may occupy is larger than for a random arrangement of spheres. At some point, the 
spheres become arrested giving rise to the divergence of the viscosity. These considerations suggest 
that the relevant quantity in expression (jlOp is not <fi but an effective filling fraction (j> e ff that 
incorporates geometrical information. We propose the following expression for 4> e ff'- 



where c is a structural or crowding factor that takes into account the arrangements of the spheres 
in the suspension and is given by 



with (p c the filling fraction at maximum packing. The main characteristic of (f>eff is that it is the 
most simple way to incorporate geometrical information having the characteristic that 4> e ff ~ <j> 
at low volume fractions and 4> e ff = 1 at maximum packing. The introduction of <j) e ff is analogous 
to the excluded volume correction proposed by van der Waals to incorporate the finite size of 
the molecules in the equation of state of the ideal gas. In fact, considering the excluded volume 
correction, the equation of state of an ideal gas is replaced by p = p e ffksT, with p e ff = pj (1 — vp). 
Here p = N/V is the number density of molecules in the volume V and v is the volume of each 
molecule. In our case, the excluded volume is due not only to the finite size of the spheres but also 
to the volume that is not accessible to other spheres due to their geometrical arrangement. 

The differential effective medium theory (DEMT) is based in a progressive addition of spheres 
to the sample. Consider a suspension of spheres with a filling fraction 4> e ff- Then, suppose that 
we increase by <5</> e // the particle concentration in the suspension of viscosity f?(0 e //) by adding 
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a small quantity A<f> e ff of few new particles. If we treat the suspension into which we add these 
particles as a homogeneous effective medium of viscosity rj((j) e ff), then the new viscosity can be 
written as 

V (4>eff + S<f>eff) = V (4>eff) U + ^^eff) , (13) 

where the added quantity of new spheres is related to Sffreff by 

^» = T^f- < 14 > 

1 - <Peff 

The normalizing factor that appears in Eq. (|14p is due to the fact that one has to remove part of 
the effective medium in order to allocate the new particles, in other words, the new fraction of 
spheres is given by 4> e ff + 6<j> e ff = <j) e ff(l — A(j) e ff) J t-A.4> e ff, from which we find A 4> e fj- Integrating 
Eq. (fT3|) we obtain 

7?(^ = % (l-0 e// )- 5/2 . (15) 



IS], Eq.®, and [19(1, Eq.©, and also to 



This expression is similar to the empirical proposals of 

other differential effective medium theories [20I]- [2^. However, our model crucially makes use of 
the geometrical information of the system from the start through </> e //. This results in an huge 
improvement of the quantitative description when compared to experimental data as will be shown 
in Section IV. Thus, Eq. (|15p is a powerful expression accounting for a variety of situations ranging 
from disordered to ordered system of spheres, or to systems with poly-dispersion of sizes. Even 
more, by an appropriate choice of the coefficient of the <p term in Eq. (|10|) it can also be applied to 
other systems ranging from emulsions to arbitrarily shaped particles. 



III. MESOSCOPIC HYDRODYNAMICS ANALYSIS OF THE MODEL 

The model of a concentrated suspension of particles formulated in the previous section essentially 
considers the superposition of dilute suspensions of particles added in a heat bath at different 
stages. At each stage, the heat bath is considered as an effective continuum medium with an 
effective viscosity depending on the volume fraction of the particles added up to the previous 
stage. The assumption that the heat bath is a continuum medium has physical implications that 
can be elucidated by performing a mesoscopic analysis of the dynamics of the system. 

In non-equilibrium conditions imposed by a stationary flow vo(r), the dynamics of a dilute "gas" 
of Brownian particles of mass m suspended in a Newtonian heat bath can be analyzed by means of 
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a Fokker-Planck equation for the single-particle probability distribution function 



31 



32 



34, 



When the volume fraction of the suspended particles is finite, hydrodynamic interactions play an 
important role and then the description is usually performed by using a N— p article distribution 



n 



function obeying a multivariate Fokker-Planck or Smoluchowski equations [15|, L 

However, to understand the physical assumptions underlying the differential model previously 
introduced, we will follow a similar procedure in which, at each stage, the new particles added con- 
stitute a dilute suspension in an effective continuum medium. Under this assumption, the dynamics 
of the system can be performed by means of the single-particle probability distribution /(r, u, t), 
depending on time t and the instantaneous position r and velocity u of the particles. Assum- 
ing that the corresponding fluctuation-dissipation theorem is valid for the actual non-equilibrium 



steady-state 



32, 



341 ] . it is possible to show that / obeys the Fokker-Planck equation 

k B TdP~ 



|p + V .(uP)=^ 



(u-v )P + 



(16) 



m du 

where V = d/dr, k B T is the thermal energy (assumed constant) and (5 is the friction coefficient 
per mass unit of a particle. Here, vo is the stationary solution of the Navier-Stokes equation 
poVoVvo = V • P° for the heat bath, with the appropriate boundary conditions on the surface 
of the particle and po is the constant density of the bath. Equations similar to (|16p have been 
derived for related systems where the external flow introduces corrections leading to the breaking 



of the fluctuation-dissipation theorem 



31 



33j, 



39l | . Here, we will assume that Eq. (p76|) is valid for 



arbitrary externally applied velocity gradients. 

Eq. (I16p allows one to calculate a set of evolution equations for the moments of the probability 



distribution, constituting the "hydrodynamic" description of the system [3ll . l37l |. Taking into 
account the usual definitions for the hydrodynamic-like fields: mass density p(r,t), velocity v(r, t) 
and pressure tensor P fc (v, t) (see reference 31]), the evolution equation for p is found to be 



dp 

at 



-V • J 



D, 



(17) 



where the diffusion current Jd(i", t) = pv obeys 38] 

1 DJ D , 



V-P A 



(18) 



Here D/Dt is the material time derivative. Eq. (|18l) contains the pressure tensor P fc of the 



suspended phase whose evolution equation is given by 



31 



37] 



2 p fc + /r 1 ^-p k 



knT 

2-2— pi - 2pD Vv s , 



(19) 
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where the symbol s stands for the symmetric part of a tensor, 1 is the unit tensor and Dq = is 
the single-particle diffusion coefficient. In writing the second term at the left-hand side of Eq. (|18p . 
we have neglected contributions of the order j3~ 2 . The last equation can be separated into two 
equations by considering P fc = pi + L, with L(r, i) the traceless stress tensor. Then, one obtains 
an evolution equation for the pressure p, and other one for the non-diagonal terms L, associated 
with the dissipation in the system. Assuming that the diagonal part relaxes faster than the other 
components, one has the state equation p = ^fp- On the other hand, the resulting evolution 
equation for L can be recasted in the form 

2L + /T^L = -2pA)(Vv )° - 2pD (VJ)°, (20) 

where we have used the definition of the local diffusion current J = Jd — p^o = p(v — vo) and 
neglected terms of the order /> _1 Vp. 

An evolution equation for the total stress tensor Q° of the suspension can be obtained from 
Eq. (|20|) by considering the symmetric traceless part of the stress tensor of the heat bath: P° = 



2//o ( Vvq) u , where rjo is the viscosity of the (Newtonian) bath [?]. Thus, multiplying Eq. (f20|) 
by —pop~~ l with p the average density of the suspended phase, adding P° and neglecting the 
convective term in the total time derivative, one arrives at the following evolution equation for the 
total traceless stress tensor of the suspension Q° = P° — 2poP~ 1 L , 

Q° + P~ l §- t Q° = 2%(W ) + 2D m ~p-\VJ)\ (21) 
where we have defined the effective viscosity coefficient r] e of the whole suspension as 

rie = ^(1 + A)/Wo~ *) > ( 22 ) 

which takes into account the effects of Brownian motion of the particles. In deriving Eqs. (|2ip and 
(|22p we have assumed that terms of the order Dqp~ 2 J ■ S/p are negligible. This is consistent with 
the approximation Ap/p <C 1 used to derive Eq. ([20]) . 

The dynamic viscosity of the suspension as a function of the frequency u and the wave vector 
k, can be calculated by linearizing the hydrodynamic equations (|17|) - (|2ip . With this purpose, we 
will assume that p = p + Sp, vq = v s + 5u, L = 5~L and Q° = Q{] + <5Q°, with 5p, 5u and 
5Q° the deviations with respect to the average values p, v s and Qq. Here, 5Q° contains the 
effect on the stresses induced by the particles when locally perturb the velocity field of the bath: 
5Q° = 2r ?0 V<5u° - 2 / 9 p- 1 <5L. 
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Then, at first order in the deviations, from Eqs. (|17|) . (|18p and (|2ip . one may derive the following 
set of fluctuating hydrodynamic equations [stJ 



^Sp + v s • V5p = - V • (S3) , 



(23) 



5 J + = -D V6p - ^P^pPo 1 ^ ■ 5Q° - 2/T x V • <5L, 



(24) 



<5Q° + /? _1 ^Q° = 2r/ e (V5u )° + 2D pp 1 (V.5J) . 



(25) 



Assuming that v s is a linear flow, then the correction 5u due to the presence of the particle is 
proportional to the constant tensor Vv s [411 ] . Following Ref. [41], Eq. (|25p can be written in the 
form 



SQ° + J^Q° = 50% (Vv s )° + 2D ppo 1 (V5J )° . 



(26) 



An equation similar to (|26p can be derived for <5L. Taking the space-time Fourier transform of 
Eqs. (|23p . (|24p and (|26p . one may solve for <5Q°. Assuming for simplicity that k = (A;, 0,0) and 
neglecting terms smaller than /3~ 2 , one obtains 

5(j)rj e 



xy 



-, ■ a-1 Do/?" 1 *: 2 ( VVs W 

1 — imp 1 — 



(27) 



Using Eq. $2} and Q° = Qg + 5Q°, we obtain a dynamic viscosity which follows a second order 
continued-fraction expansion 



r](u, k; <f)) ~ r/o 



i + 5 



1 + So' 1 



1 LLU P (l-iuj/3- 1 ) 



(28) 



where we have introduced the Schmidt number Sc = %(-DoPo) _1 ) which for typical suspensions is 
usually a small quantity < 10~ 3 . By expanding Eq. ([28]) up to first order in u and second order 
in k, we obtain an expression for the viscosity similar to that reported in Ref. [7[. This (w, /in- 
dependent correction could be significative for very small times t ~ Therefore, at (diffusion) 
times t 3> Eq. (|28j) reduces to Einstein's formula: rj(<j)) ~ % [l + §^>1 . 
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A. Generalization to finite volume fractions 

To generalize results obtained in the previous section to the case of finite volume fractions, 
we will follow arguments similar to those of Section II, but applied to the long-time limit of the 
evolution equations derived in the previous section. Then, the subindex % will be used to indicate 
the stage % in which, according to Einstein's relation, the effective viscosity is rj((pi) ~ 770 [l + §0j] , 
with the volume fraction at stage i, (pi, defined as in Eq. (jlll) . 

Thus, assuming that the suspension of particles at stage i plays the role of an effective (contin- 
uum) medium, it is clear that the probability distribution of the new phase (i + 1) will satisfy a 
Fokker-Planck equation similar to Eq. (|16p . where now the friction coefficient can be denoted by 
and the diffusion coefficient by Dq +1 ^ = kgT '/ '(m/3j + i). Following a similar procedure as the 
already indicated in this section, from the Fokker-Planck equation at stage i + 1 the linearized set 
of evolution equations for the fluctuations of the hydrodynamic fields is 

d 

— 5p i+1 + v s • V5 Pl+1 = -V • (tfJi+i) , (29) 
5J l+1 + fc+iQ- t Mi+i = -Do +1 V5p l+1 - -p- + \p i+1 p^V ■ 6Q° +1 - 2/3-\V • «JL i+1 , (30) 

<5Q°+i + /T+i J^Q°+i = 2 Ve>i (V5u)° + 2Dl +1 ~ Pl+1 p G 1 (VS3 i+1 )° . (31) 

where we have used the definition Q° +1 = Q° — 2p p~^ 1 'Li + i. In writing Eq. (|3"T|) an important 
assumption has been done: We have assumed that the traceless stress tensor Q-* of the effective 
medium has relaxed to its stationary value, thus implying that only the relaxation time (3^ 1 enters 
into Eq. (I3ip . This hypothesis allows us to obtain the following relation for the effective viscosity 
of the system at the stage i + 1: 

(32) 

If we now take into account that the total volume fraction occupied by the particles at the stage 
i + 1 is given by <p = (f>i(l — fa+i) + <j>i+i an d assume that at each stage a very small (differential) 
quantity of particles is added, then: cpi + i = 6(j)/(l — (pi), with 5<f> = (j> — (pi. Substituting this 
relation into Eq. (|32p . we finally obtain the differential relation Srj/rji ~ (5/2)50 / (1 — (pi) with 
5rj = rji + i — rji that after integration leads to 

r? (</>)= r/o (l-</>r 5/2 , (33) 



i+l, 



Vi{(Pi 



1 + 
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which has the same form of Eq. (|15p . As in Section II, it must be emphasized that geometrical 
information must be introduced by taking into account that Eq. ([33]) depends upon the effective 
volume fraction cj> e ff defined in Eq. (fTTI) . 



IV. RESULTS AND COMPARISON WITH EXPERIMENTS 



In this section we compare the predictions of the improved DEMT with several experimental 
results and with the results obtained from other theories. To carry out these comparisons, it is 
important to notice that the value of 4> c in Eq. (|12p is a free parameter of the theory to be chosen in 
order to best fit the experimental results. Nonetheless, this parameter can be chosen beforehand 
based on physical arguments, and then used to compare with specific experimental situations. 

This is done in Fig. 1, where the behavior of the relative viscosity r\ (4>) /rjo is compared with 



experimental results of de Kruif et al. 



43j and of Krieger 



451 ] for high and low-shear rates as a 
function of the volume fraction <p. The lower curve is the prediction of Eq. (|15p with <p c = 0.7404, 
which corresponds to FCC close packing. This gives an excellent agreement with the experimental 
results for the case of high-shear rates. The upper curve represents the prediction of Eq. (|15H 
with (p c = 0.63, which corresponds to the random close packing (RCP) of identical spheres. The 
comparison with the experimental results for low-shear rates is again excellent. These results are 
consistent with the known fact that for low-shear rates the spheres remain disordered while at high- 
shear rates the equilibrium microstructure of the dispersion is completely destroyed and the spheres 
adopt an ordered FCC configuration. In Fig. 2 we use again the value 
theory with the experimental results obtained by Saunders 
by Kops-Werkhoven and Fijnaut 



= 0.63 to compare our 
44i | , by Krieger and coworkers 4f]] , and 
46p . We also plot the calculated values obtained by Beenakker[jJ] 



with a theory that fully takes into account the many-body hydrodynamic interactions between an 
arbitrary number of spheres. Notice again that the agreement with the experimental results is 
excellent over the entire fluid range. 

In Fig. 3 we compare r? (<j)) /r]o with the reduced viscosity measurements by van der Werff et 
al. [12I and by Zhu et al. 



and Felderhof 



481 1 at high frequencies. We also plot the values obtained by Cichocki 

1 n 

30l | and the predictions using a proposal by Verberg et al. [15] which considers 



the fraction of colloidal particle pairs at contact, Eq. (|7|). We used as a parameter to fit the data 
(j) c = 0.8678. We note that Verberg's predictions lie below the experimental points while the results 
of the improved DEMT fits the data very well. Similarly, Fig. 4 shows the zero frequency low-shear 
viscosity as obtained by our model using <j> c = 0.63, and it is compared with experimental results 
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and the model of Verberg et al. [151 ] . Again, a very good agreement is found. 

A more sensitive representation of the relative viscosity data is the one proposed by Bedeaux 



3 



0(1 + 5(0)), (34) 



V (<t>) A?0 + 2 

where S(<p) is an unknown function of the volume fraction. In Fig. 5, we represent the function 
S((f)) for various experimental results obtained in a variety of experimental situations and compare 
with the values obtained with our model. We take the same values of 4> c used previously for the 
low- and high-shear cases. The agreement with the experimental data is again very good. The 
last set of experimental points corresponds to the high-frequency limit of the real part of the 
complex viscosity as obtained by the linear dynamic measurements performed by van der Werff 



and coworkers 



12]]. The best fit to the real part, obtained with C = 0.8678 gives an excellent 



agreement with the experimental values. As a comparative, we also show the predictions for the 
righ- frequency limit, obtained by Beenakker [9], by Batchelor and Green J], by Russel and Gast 
lGj, and by Verberg et al. [15]. In this very sensitive representation it is seen that only the 
improved DEMT matches the data at all volume fractions. 

Finally, let us notice that the functional form for the reduced viscosity in the improved DEMT 
as given by Eq. (|15p makes no reference to any material parameter neither to the geometrical 
information of the system. This means that a plot of the reduced viscosity i](<j)) /r/o vs. e // 
should be independent of the experimental details or of the sample. This is confirmed in Fig. 6 
where we plot [rj (0) /rjo]~ 2 ^ vs - 4>eff f° r & h the data sets considered in the previous figures. The 
straight line is the prediction of the improved DEMT, Eq. (|15p . As can be seen, the generality of 
the correlation and the agreement with our model is quite remarkable. 

Now we proceed to compare our expression for the viscosity with the functional form of other 
theories. First, let us notice that putting (ft = <ft c in the denominator of Eq. ljlip and replacing the 
'actor 5/2 by 2 in Eq. (|10p we obtain the well known empirical expression obtained by Quemada 
Jjj], Eq.(|9]). This expression although gives the right divergence at high filling fractions does not 
have the correct low volume fraction limit given by Einstein's expression. 

On the other hand, putting 4> e ff = 0/0c in the denominator of Eq. (|14p we recover the empirical 



expression obtained by Krieger and Dougherty 



Eq.(|HJ). Nonetheless the virial expansion of this 



22[ | . A generalization 



expression gives values for the second virial coefficient, &2, that are too small 
of this model but including a dependence of the crowding factor on the volume fraction was proposed 
by van de Ven 22| • However this model has the shortcoming that gives rise to an expression for the 
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viscosity that is more involved and it includes a second adjustable parameter besides (f) c which is 
chosen in order to adjust the second virial coefficient to the values predicted by other theories. In 
contrast, our proposal does not contain adjustable parameters besides 4> c and the series expansion 
gives the right values for the second virial coefficient in a natural way as it is shown below. 

Expanding Eq. (ll5p in power series of (j> for specific choices of <f> c allows one to compare the 
improved DEMT with other theories. Choosing the RCP value for the volume fraction at maximum 
packing, <j) c = 0.63, we obtain up to O(0 3 ) 



which compares very well up to second order with the result that Cichoki and Felderhof [30] have 
obtained on the basis of the pair Smoluchowski equation which are exact 0(<f> 2 ). With Brownian 
motion contributions they find for the second virial coefficient k<i = 5.91 [see Eq.([5])]. Batchelor 
and Green [J] including Brownian motion have obtained k<i = 6.17, as shown in Eq.([2]). Our result 
differs form the exact one by 1.2%, while the result by Batchelor and Green differs by 4.4%. 
Taking FCC close packing cp c = 0.7404, and expanding Eq. (fT5|) up to 0((j) 3 ) we obtain 



which compares well up to second order with the result by Bedeaux and collaborators [7j], &2 = 4.8, 
as shown in Eq.([3]) and even better with the result of Cichocki and Felderhof without including 
Brownian motion k2 = 5.00, as shown in Eq.([4|). Our result differs form the exact one by 5%, while 
the result by Bedeaux by 4%. 

Using the value that best fits the high-frequency regime <p c = 0.8678, one finds for the real part 
of the complex relative viscosity 



which compares well with the result by Batchelor and Green [4|, ki = 5.2, as shown in Eq.([T]) and 
in which Brownian motion is not taken into consideration. 

Finally, it is interesting to find a value for cj) c that best fits the theoretical results obtained by 
Verberg et al. [ij]]. This value is cp c = 0.95238 and produces the following series expansion up to 
O(cf) 7 ), just to show a few terms 






rj (</>) /rj = 1 + ^<f> + 4.5<f> 2 + 7.OO60 3 + 1O.O410 4 + 13.63580 5 + 17.8297/ + 22.66730 7 + 0(/), 
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which is exact up to 0{4> 2 ) with the series expansion of the result by Verberg et al. which reads 

rj (<f>) /rj = 1 + ^<f> + 4.5^ 2 + 7<j) 3 + 10^ 4 + 13.50 5 + 17.50 6 + 220 7 + 0(<fP), 
and the higher order terms are also amazingly close. 

V. CONCLUSIONS 

In summary, we have formulated an improved differential effective medium theory that accounts 
for the dependence of the viscosity of hard-spheres suspensions at arbitrary volume fractions. This 
theory takes into consideration geometrical information of the system from the start through an 
effective filling fraction and assumes a hierarchy of relaxation times, as we have shown by performing 
a mesoscopic analysis based on Fokker-Planck equations. 

The agreement with experimental data at low- and high-shear rates and in the high-frequency 
limit is remarkable and superior to other theories. Also, the virial expansion of our model compares 
well with the second virial coefficient obtained by a number of different theories in all the regimes 
considered. In the particular case of the high-frequency limit, it is possible to find a value for <p c 
for which the virial expansion is amazingly close to higher orders to that of the model by Verberg 
et al. Although our proposal can be generalized to describe the viscoelastic behavior at lo / 
by simply allowing <j) c to take complex values, we have not a derivation based on simple physical 
arguments to justify this generalization. Such a goal requires further research. 
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FIG. 1: Relative viscosity rj (<f>) / 770 at low and high shear rates as a function of the volume fraction <fi. 
The lines correspond to the predictions of the improved DEMT with RCP <f> c = 0.63 (upper line) and 
close packing at FCC (f> c — 0.7404 (lower line). The measured data are from Refs. [43] (circles) and 
(triangles). 
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FIG. 2: Volume fraction dependence of the reciprocal of the relative viscosity rjo/r] (<£). The solid line is the 
result of the improved DEMT with <p c = 0.63. The dashed line is taken from Ref. [9(. The measured data 
are from refs. [44| (squares), 45 1 (triangles), and [46| (circles). 
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FIG. 3: Relative infinite-frequency viscosity r\ (cj>) / 770 as a function of the volume fraction cf>. Squares, 
Zhu et al. circles, van der Werff et al. triangles. Cichocki and Felderhof {3^. The dashed line 



corresponds to the prediction by Verberg et al. 
with <h r = 0.8678. 



15j , and the solid line is the result of the improved DEMT 
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FIG. 4: Relative zero-frequency viscosity rj ((f)) / 770 at low shear rate as a function of the volume fraction 
4>. Triangles, van der Werff et al. Q; crosses, van der Werff and de Kruif circles, Jones et al. 
squares, Papir and Krieger [5o| . The dashed line corresponds to the prediction by Verberg et al. and 
the solid line is the result of the improved DEMT with <fi c = 0.63. 



21 




FIG. 5: Representation of various viscosity data as suggested by Bedeaux (Ref. [471]). Squares, SJ18 low- 
shear limit (Refs. H and Q); 

circles, SJ18 high-shear limit (Refs. |13j and [3]); triangles, high-frequency 
limit of the real part of the complex shear viscosity (Ref. 12j). Thick solid lines are the results of the 
improved DEMT with <\> c — 0.63 (upper), <fi c = 0.7404 (medium), and 4> c — 0.8678 (lower); thin dashed line 
is Beenakker's result (Ref. Q), thin dotted line is Verberg et al. prediction (Ref. [la]), thin dash-dot line is 
the result of Russel and Gast (Ref. [l0(), and the thin dash-dot-dot line is the prediction of Batchelor and 



Green (Ref. £|) 
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FIG. 6: Relative viscosity [77 (<f>) /770] 2 ^ 5 as a function of the effective volume fraction cj> e ff- Symbols 
correspond to all the data sets of the previous figures and the solid line corresponds to the prediction of the 
improved DEMT, Eq.fH]). 



